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| Abstract. We present an accurate methodology for representing the physics of waves, 

Y\ for periodic structures, through effective properties for a replacement bulk medium: 

• ^h This is valid even for media with zero frequency stop-bands and where high frequency 

O | phenomena dominate. 

Since the work of Lord Rayleigh in 1892, low frequency (or quasi-static) behaviour 

c/3 has been neatly encapsulated in effective anisotropic media; the various parameters 

. — i come from asymptotic analysis relying upon the ratio of the array pitch to the 

wavelength being sufficiently small. However such classical homogenization theories 
r^ break down in the high-frequency, or stop band, regime whereby the wavelength to 

Oh pitch ratio is of order one. Furthermore, arrays of inclusions with Dirichlet data lead to 

a zero frequency stop band, with the salient consequence that classical homogenization 
^~H is invalid. 

J^ Higher frequency phenomena are of significant importance in photonics (transverse 

qq magnetic waves propagating in infinite conducting parallel fibers), phononics (anti- 

£" — plane shear waves propagating in isotropic elastic materials with inclusions), and 

**| platonics (flexural waves propagating in thin-elastic plates with holes). Fortunately, 

the recently proposed high-frequency homogenization (HFH) theory is only constrained 

by the knowledge of standing waves in order to asymptotically reconstruct dispersion 

curves and associated Floquct-Bloch eigenfields: It is capable of accurately 

L* representing zero- frequency stop band structures. The homogenized equations are 

. ^h partial differential equations with a dispersive anisotropic homogenized tensor that 

/\ characterizes the effective medium. 

3 We apply HFH to metamaterials, exploiting the subtle features of Bloch dispersion 

curves such as Dirac-like cones, as well as zero and negative group velocity near stop 
bands in order to achieve exciting physical phenomena such as cloaking, lensing and 
endoscope effects. These are simulated numerically using finite elements and compared 
to predictions from HFH. An extension of HFH to periodic supercclls enabling complete 
reconstruction of dispersion curves through an unfolding technique is also introduced. 

PACS numbers: 41.20.Jb,42.25.Bs,42.70.Qs,43.20.Bi,43.25.Gf 



1. Introduction 

Over the past 25 years, many significant advances have created a deep understanding 
of the optical properties of photonic crystals (PCs) (Yablonovitch 1987, John 1987); 
such periodic structures prohibit the propagation of light, or allow it only in certain 
directions at certain frequencies, or localize light in specified areas. This sort of 
metamaterial (using the consensual terminology for artificial materials engineered to 
have desired properties that may not be found in nature, such as negative refraction, see 
e.g. (Smith et al. 2004, Ramakrishna 2005)) enables a marked enhancement of control 
over light propagation; this arises from, for instance, the periodic patterning of small 
metallic inclusions embedded within a dielectric matrix (Pendry et al. 1996, Nicorovici 
et al. 1995). PCs are periodic devices whose spectrum is characterized by photonic 
band gaps and pass bands, just as electronic band gaps exist in semiconductors: In 
PCs, light propagation is disallowed for certain frequencies in certain directions. This 
effect is well known (Joannopoulos et al. 2008, Zolla et al. 2005) and forms the basis 
of many devices, including Bragg mirrors, dielectric Fabry-Perot filters, and distributed 
feedback lasers; all of these contain low-loss dielectrics periodic in one dimension, so 
are one-dimensional PCs. Such mirrors are tremendously useful, but their reflecting 
properties critically depend upon the frequency of the incident wave in regard with 
its incidence (Markos & Soukoulis 2008). For broad frequency ranges, one wishes 
to reflect light of any polarization at any angle (which requires a complete photonic 
band gap) and for Dirichlet media (i.e. those composed with microstructure where the 
field is zero on the microparticles) such a gap occurs at zero-frequency It is therefore 
possible to create seismic shields for low-frequency waves (Brule et al. 2013, Antonakakis 
et al. 2013b). Extending these ideas to higher dimensions allows for a much greater range 
of optical effects: all-angle negative refraction (Luo et al. 2002, Zengerle 1987), ultra- 
refraction (Farhat et al. 2010, Craster et al. 2011) and cloaking at Dirac-like cones (Chan 
et al. 2012): We will demonstrate here that an effective medium can be created that 
reproduces these effects. 

The considerable recent activity in this area is, in part, fuelled by advances 
in numerical techniques, based on Fourier expansions in the vector electromagnetic 
Maxwell equations (Joannopoulos et al. 2008), Finite Elements (Zolla et al. 2005) and 
Multipole and lattice sums (Movchan et al. 2002) for cylinders, to name but a few, that 
allow one to visualize the various effects and design PCs. The multipole method takes 
its root in a seminal paper of Lord Rayleigh (Rayleigh 1892) which is the foundation 
stone upon which the current edifice of homogenization theories is built. More precisely, 
in that paper, John William Strutt, the third Lord Rayleigh, solved Laplace's equation 
in two dimensions for rectangular arrays of cylinders, and in three dimensions for cubic 
lattices of spheres. However, a limitation of Rayleigh's algorithm is that it does not hold 
when the volume fraction of inclusions increases. Multipole methods, in conjunction 
with lattice sums, overcome such obstacles and lead to the Rayleigh system which is 
an infinite linear algebraic system; this formulation, in terms of an eigenvalue problem, 
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Figure 1. High Frequency Homogenization (HFH) principle: An elementary cell (a) of 
sidelcngth 21, modelled by a fast oscillating variable £, is repeated periodically within 
a supercell (b) of sidelcngth L, modelled by a slow variable X, which is itself repeated 
periodically in space (c). One then assumes that the parameter e = l/L is small, and its 
vanishing limit thereafter studied. The leading order, homogenized, term of Floquet- 
Bloch eigenfields within the crystal are then sought as Uo(X, £) = /o(X)[/o(£; ^o)> 
wherein /o accounts for variations of the fields on the order of the supercells, and Uo 
captures their fast oscillations in the much smaller cells, when either periodic or anti- 
periodic conditions are enforced on the cells: Perturbing away from these standing 
waves of frequency Qq allows for a complete reconstruction of the Bloch spectrum and 
associated Floquet-Bloch eigenfields. 



facilitates the construction of dispersion curves and the study of both photonic and 
phononic band-gap structures. In the limit of small inclusion radii, when propagating 
modes are very close to plane waves, one can truncate the Rayleigh system ignoring 
the effect of higher multipoles, to produce a series of approximations each successively 
more accurate to higher values of filling fraction. At the dipole order, one is able to fit 
the acoustic band, which has a linear behaviour in the neighbourhood of zero frequency, 
except of course in the singular case whereby Dirichlet data is enforced on the boundary 
inclusions (Poulton et al. 2001) and there is a zero frequency stop-band: This is a 
substantial limitation of this truncation. We will derive here the exact solution for 
zero radius Dirichlet holes, avoiding multipole expansions, and show dispersion curves 
whose features lead to interesting topical effects, such as Dirac-like cone cloaking and 
degenerate points (Chan et al. 2012) which emerge very naturally in this exact solution. 
Although the numerical approaches discussed briefly are efficient they still do 



require substantial computational effort and can obscure physical understanding and 
interpretation. Here we substantially reduce the numerical complexity of the wave 
problem using asymptotic analysis; this has been developed over the past 35 years 
by applied mathematicians primarily for solving partial differential equations, with 
rapidly oscillating periodic coefficients, in the context of thermostatics, continuum 
mechanics or electrostatics (Bensoussan et al. 1978, Jikov et al. 1994, Milton 2002). 
The available literature on such effective medium theories is vast, but it seems that 
only a very few groups have addressed such problems as the homogenization of media, 
with moderate contrast in the material properties, in three-dimensions (Guenneau & 
Zolla 2000, Wellander & Kristensson 2003) and high-contrast two-dimensional (Zhikov 
2000, Bouchitte & Felbacq 2004, Cherednichenko et al. 2006) photonic crystals, that 
have important potential applications in photonics. Besides this, the aforementioned 
literature does not address the challenging problem of homogenization for moderate 
contrast photonic crystals near stop band frequencies: Classical homogenization is 
constrained to low frequencies and long waves in moderate contrast PCs (Silveirinha & 
Fernandes 2005) and so-called high-contrast homogenization only captures the essence 
of stop bands in PCs when the permittivity inside the inclusions is much higher than 
that of the surrounding matrix (the contrast being typically on the order to e -2 where 
e is the array pitch which in turn is much smaller than the wavelength). This latter 
area of homogenization theory is fuelled by interest in artificial magnetism, initiated by 
the work of O'Brien and Pendry (O'Brien & Pendry 2002). However, moderate contrast 
one-dimensional photonic crystals have been recently shown to display not only artificial 
magnetism, but also chirality (Liu et al. 2013), which is an onset to negative refraction 
(Pendry 2004) . 

For all these reasons, there is a strong demand for high-frequency homogenization 
theories of PCs in order to grasp, and fully exploit, the rich behaviour of photonic 
band gap structures (Notomi 2000, Paul et al. 2011). This has created a suite 
of extended homogenization theories for periodic media called Bloch homogenisation 
(Conca et al. 1995, Allaire & Piatnitski 2005, Birman & Suslina 2006, Hoefer & 
Weinstein 2011). There is also a flourishing literature on developing homogenized elastic 
media, with frequency dependent effective parameters, also based upon periodic media 
(Nemat-Nasser et al. 2011): There is considerable interest in creating effective continuum 
models of microstructured media that break free from the conventional low frequency 
homogenisation limitations. 

In this paper, we apply the HFH theory proposed by one of us three years ago 
(Craster et al. 2010) to the physics of zero frequency band gap structures, not only in 
the context of photonics, but also phononics for anti-plane shear waves in periodic arrays 
of inclusions, and platonics with flexural waves in pinned plates. The latter analysis is 
made possible thanks to the extension of HFH to plate theory, developed by two of us 
two years ago (Antonakakis & Craster 2012). Our aim is to show the universal features 
of stop band structures thanks to HFH, and to further exemplify their potential use in 
control of light and mechanical waves, with novel applications ranging from cloaking 




Figure 2. Dynamic effective anisotropy: A time-harmonic source located inside a 
square array of pitch 2 consisting of 196 circular constrained holes of radius 0.4 leads 
to a wave pattern resembling a Saint- Andrews' cross at frequency ft = 1.966 (a) FEM 
(b) HFH, and a Saint-George's cross at frequency = 2.75 (c) FEM (d) HFH. 



(Milton & Nicorovici 2006, Guenneau et al. 2007) to anti-earthquake seismic shields 
(Antonakakis et al. 2013a, Brule et al. 2013). 

We show in Fig. [T] what HFH does in practice: It focuses its attention on the 
physics within a supercell of sidelength L, which is the long- scale, and further captures 
the fine features of field's oscillations inside an elementary cell of sidelength 2/ much 
smaller than L (Craster et al. 2010). The wavelength need not be large compared to 
/ in order to perform the asymptotic analysis, as the small parameter e = l/L only 
requires the supercell to be much larger than its constituent elementary cells; this 
contrasts with classical homogenization that assumes the cells are much smaller than 
the wavelength. In this way, one replaces a periodic structure by a homogenized one, on 
the long-scale, at any frequency and classical homogenization is just a particular case of 
HFH (Antonakakis et al. 2013a). An illustration of the HFH theory is in Fig. [2] the left 
panels are from full finite element simulations of an array of cylindrical holes (infinite 
conducting boundary conditions of the holes which are for transverse magnetic (TM) 
waves in electromagnetism, or clamped shear horizontal waves in elasticity) excited at 
the array centre by a line source and the right panels replace the array with its effective 
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Figure 3. Panel (a) shows a square array of circular cylinders, with the elementary 
cell as the dashed box. In (b) the irreducible Brillouin zone is shown together with 
lettering necessary for our discussion of folding in section |4.1.2| 



HFH medium which reproduces the large scale behaviour. For the frequencies chosen, 
two very distinct cross-shaped propagations, a Saint Andrews' cross (or saltire) Fig. 
|2^a),(b), and a Saint George's cross Fig. [2](c),(d) are achieved; these distinctive shapes 
are created by excitations at precise frequencies that are identified from the dispersion 
curves, (Craster et al. 2012), this strong frequency dependent anisotropy is a recurrent 
feature of PCs. One sees that the HFH neatly reproduces in Fig. |2](b,d) the fine features 
of the full finite element simulations. Interpretation is provided by using the Brillouin 
zone of Fig. [3] and dispersion curves shown in Fig. |4j The square array of holes behaves 
as two dramatically contrasting effective media at frequency Vl = 1.966, see Fig. p[a,b), 
which stands on the edge of the zero-frequency stop band that coincides with the lower 
edge of the second stop band; and at frequency Q = 2.75 which stands on the upper 
edge of the second stop band, wherein the dispersion curve is nearly flat is the MX 
direction; both these frequencies are well beyond the range of applicability of classical 
homogenizat ion . 

Our aim is to generate HFH for arrays of holes with Dirichlet conditions, so this is for 
TM waves in electromagnetism, and then compare through full numerical simulations 
of HFH and finite elements to show how the behaviour of the numerous interesting 
optical effects emerge naturally through the coefficients of our effective medium. This 
is intertwined with a knowledge and understanding of the wave structure for a perfect 
medium where the dispersion relations connecting the phase shift across an elementary 
cell with the frequency is key. In section [2] we set up the mathematical framework 



leading to the effective medium equation (15), this is in a non-dimensional setting so for 



clarity in section [3] the main results are summarized in a dimensional setting. There are 
degenerate cases that are of substantial interest and connect with Dirac-like dispersion 
(Chan et al. 2012) and these lead to a different effective equation also outlined in section 

EJ 



A first step in verifying the HFH theory is by asymptotically reproducing the 
dispersion curves for perfect infinite arrays, and we use both circular and square holes 
for illustration, as shown in section |4j Local to the edges of the Brillouin zone there 
are standing wave frequencies and the local asymptotic behaviour is given from the 
effective medium equation. Importantly, one can then predict apriori, from the signs 
of the HFH coefficients with the effective equation, how the bulk medium will behave. 



As a further illustration, holes that degenerate to a point are treated (section 4.3) as 
there is then an exact and simple solution for the dispersion relation, furthermore in 
this case degeneracies occur. Finally, in section [5] we compare HFH with full numerical 
simulations showing the full range of features available and how HFH represents them. 
Some concluding remarks are drawn together in section |6j 

2. General theory 

A two dimensional structure composed of a doubly periodic, i.e. periodic in both x and 
y directions, square array of cells with, not necessarily circular, identical holes inside 
them is considered (see Fig. [3]). We are primarily concerned with electromagnetism 
where Maxwell's equations separate into transverse electric (TE) and TM polarizations 
and the natural boundary conditions on the holes are then Neumann and Dirichlet 
respectively; the asymptotics for the former are considered in (Antonakakis et al. 2013 a) 
and applications in metamaterials are illustrated with split ring resonators (Pendry 
et al. 1999). In the current article our emphasis is rather different, and is upon TM 
polarization for which the Dirichlet conditions induce different effects, such as zero- 
frequency stop bands, and require a modified theory. 

Assuming an exp(— iut) time dependence that is considered understood, and 
henceforth suppressed, the governing equation is, 

V x • [o(£)V x w(x)] + w 2 p(£Mx) = (1) 

on — oo < x±, x 2 < oo where to is the angular frequency. In the periodic setting, each cell 
is identical and the material is characterized by two periodic functions on the short-scale, 
in £ = (xi/l,X2/l), namely o(£) and p(£). In the context of photonics, these are the 
inverse of permeability (p _1 ), and permittivity (e), respectively, which are related to the 
wavespeed of light in vacuum c via the refractive index n as follows: efi = n 2 c~ 2 . The 
unknown u physically is the longitudinal component E z of the out-of-plane electric field 
E = (0,0, E z ), bearing in mind that the transverse magnetic field H = (H x ,H y ,0) is 
retrieved from Maxwell's equation H = iw _1 p _1 V x E. In the context of phononics, a(£) 
and p(£) would stand respectively for the shear modulus and the density of an isotropic 
elastic medium, and the unknown u would correspond to the out-of-plane displacement 
(amplitude of SH shear waves). One can also easily draw analogies with acoustic pressure 
waves in a fluid, or linear surface water waves, which explains the activity related to 
finding correspondences between models of electromagnetic (Ramakrishna 2005) and 
acoustic (Craster & Guenneau 2013) metamaterials. The length of the direct lattice 



8 

base vectors are taken equal and to be 21, and define a short lengthscale. The overall 
dimension of the structure is of a much longer length-scale, L; the ratio of scales, 
e = l/L <C 1 then provides a small parameter for use later in the asymptotic scheme. 

A non-dimensionalization of the physical functions by setting, a = aoa($,) and 
p = pop($,) is convenient. Equation ([I]) then becomes 

^ 2 V X • [a(£)V x tt(x)] + tt 2 p(£)u(x) = with n = ^- (2) 

Co 



as the non-dimensional frequency and Co = v^o/Po- The two scales, /, L, then motivate 
two new coordinates namely X = x/L, and £ = x/7 that are treated as independent, 
and placing this change of coordinates into (pi) we obtain, 

V r [a(|)V^(X,|)] + n 2 p(0«(X^) 

+ e[2o(£) V € + V € a(£)] ■ V X «(X, £) + e 2 a(£) V^u(X, £) = 0. (3) 

For wave propagation through a perfect lattice there exist standing waves at specific 
eigenfrequencies, f2 , and their eigenmodes satisfy periodic, or anti-periodic (out-of- 
phase) boundary conditions on the edges of the cell: We obtain asymptotic solutions 
based upon these standing waves. 

A key point, from this separation of scales, is that the boundary conditions on the 
short-scale are known. For instance, periodic conditions in £ on the edges of the cell are 

u|&=i = w|&=-i and u,&|&=i = u,&|&=-i (4) 

where u^ r denotes partial differentiation with respect to £j. An almost identical analysis 
holds near the anti-periodic standing wave frequencies: We illustrate the periodic case 
for definiteness. On the long-scale we set no boundary conditions, and indeed the 
methodology we present can be used even when the problem is not perfectly periodic 
(Craster et al. 2011). 

The following ansatz is taken 

u{X, £) = u (X, £)+«*i(X, £)+e 2 u 2 (X, £)+. . . , Q 2 = fig+efi?+e 2 fil+. . . :(5) 

Since w(X, £) is periodic in £ so are the Wj(X, £)'s. This leads to a hierarchy of equations 
in increasing powers of e: The leading order followed by the first order, second order 
and so on. The first three equations in ascending order yield 

( au o&)& +tt 2 pu = (6) 

( au Ui)& + Q oP u i = -( 2aM o,6 + a^u ),Xi ~ tt 2 lP uo (7) 

(««2,&)a +^lpU2 

= -auo,XiXi - (2aui,& + a&ui),*. - Vt\pui - Vl 2 2 puQ. (8) 

We start by expressing a solution for the leading order equation. For a specific eigenvalue 
flo there is a corresponding eigenmode Uo(£;flo), periodic in £, and the leading order 
solution is 

uo(X,£)=/ (X)L/o(6fto). (9) 
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The function /o(X) is determined later on, indeed the whole point is to deduce a long- 
scale continuum partial differential equation for f entirely upon the long-scale. We 
begin with a treatment of isolated eigenvalues, however repeated eigenvalues also arise, 



and require modifications, and will be dealt with in sections 2J2 and 2^ These are 
relevant to Dirac-like cones (Chan et al. 2012). 

Dirichlet conditions on the inside boundary of the cell, dS2, where S denotes the 
surface of the cell in the £ coordinates, yield 

u(X,£)| a52 = «=* «i(X,£)|as a =0, i e N (10) 

and are set in the short-scale $, so, for i = 0, Uo($; Qo)\ds 2 = 0- 

To deduce the long-scale PDE for /o(X) we follow (Craster et al. 2010), making 
changes where necessary due to the change in boundary conditions. Equation fifty is 
multiplied by U and integrated over the cell's surface to obtain 

/ / (U (aui > z.)& + SllpUoUi) dS 

= -fox, J J '{aVD&dS - ml J J P U*dS. (11) 

The first integral of the right hand side is converted to a path integral along dS 
using a corollary of the divergence theorem. The periodic conditions on dS\ and 
the homogeneous Dirichlet conditions of Uq on dS2 make it vanish. We continue by 
subtracting the integral over the cell of the product of equation (|6]) with Ui/fo to obtain 

{U {au Ui )^- Ul {aU^)^dS = -m\ J Jp(£)U*dS. (12) 




's 
Using Green's theorem equation (fT2J) becomes 



L°®{ u °^-^) ds - f ° Q2i JL mu ° dS 



(13) 



where OS = dSi + dS 2 . Due to periodic (or anti-periodic) boundary conditions of u 
and its first derivatives with respect to £j on dSi, and due to the homogeneous Dirichlet 



boundary conditions of the Wj's on dS2, the left hand side of equation (13) vanishes. It 
follows that Qi = 0, which is an important deduction as it implies that the asymptotic 
behaviour of dispersion curves near isolated eigenfrequencies is at least quadratic. We 
then solve for «i(X, £) and obtain 

ui(X,£) = /i(X)C/ (£; fto) + V x /o(X) • Uj(0- (14) 

The first term is simply a homogeneous solution that is irrelevant, and the vector 
function Ui is a particular solution. Quite often exact solutions for u\ or Ui are not 
possible to obtain so one must use numerical solutions and this is described in section 
5J We now move to the second order equation to compute /o(X) and fl|- Similarly 
to the first order equation, we multiply equation ([8| by U$ then subtract the product 
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of equation (J7|) with uij /o, finally we integrate over the cell's surface to obtain the 
following effective equation for f 

T i jfo,x t x J +ttlf = 0, with Tij = ( . ( . tJ for i,j = 1,2, (15) 

J JsP u o ab 

which is the main result of this article. In particular, the coefficients T^ encode the 

anisotropy at a specific frequency and the tjj's are integrals over the small-scale cell 

tu = [ [ aUldS + 2 f [ aU UA U dS + [ [ a^U u U dS, (16) 

J Js J Js J Js 

Uj = 2J j aU liAi U dS + J j a&V h U Q dS for i^j, (17) 

where U\ i is the ith component of vector function Ui. Note that there is no summation 



over repeated indexes for tu. The PDE for / , equation (15), is crucial as the local 
microstructure is completely encapsulated within the tensor T^; these are, for a specific 
structure at an Qq, just numerical values as illustrated in table [2} Notably the tensor 
can have negative values, or components, and it allows one to interpret and, even more 
importantly, predict changes in behaviour or when specific effects occur. The structure 
of the tensor depends upon the boundary conditions of the holes, the results here, for 
instance, are different from those of the Neumann case (Antonakakis et al. 2013a). 
Another key point is that numerically the short scale is no longer present and the PDE 



(15) is simple and quick to solve numerically, or even by hand. 

One primary aim of the present article is to deduce formulae for the local behaviour 
of dispersion curves near standing wave frequencies as these shed light upon the physical 
effects observed. For this one returns to the perfect lattice and then Floquet Bloch 
boundary conditions on the elementary cell lead immediately to /o(X) = exp(iKjXj/e). 
In this notation Kj = Kj — dj and dj = 0, tt/2, — 7r/2 depending on the location in the 



Brillouin zone. Equation (15) gives 



Q ~ Q o + tJt^j (18) 

and these locally quadratic dispersion curves are completely described by Q Q and the 
tensor Ty. 

2.1. Classical singularly perturbed zero-frequency limit 

It is natural, at this point, to contrast with usual homogenization theories. As is well- 
known (Nicorovici et al. 1995, McPhedran & Nicorovici 1997) one cannot homogenize 
the TM polarized case because the conventional approach only works at low frequencies 
in a quasi-static limit. Some authors attribute this to a singular perturbation, or non- 
commuting limit, problem whereby letting first the Bloch wavenumber and then the 
frequency tend to zero, or doing it the way around, leads to a different result (Poulton 
et al. 2001). The approach used in this article overcomes this by releasing the theory 
from the low-frequency constraint. For TE polarization (Antonakakis et al. 2013a) one 
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can explicitly connect the theories and show that the low frequency theory is a sub- 
set of the high frequency approach. The TM case differs fundamentally as there is a 
zero-frequency stop-band. 

We now prove that the usual homogenization cannot give the asymptotic behaviour 
even at low-frequency Setting Q 2 = t 2 VL 2 + . . ., that is, going to low frequency, 
and assuming that Uo(X, £) = /(X)[/o(£) is non-zero one arrives at a contradiction: 
At leading order Uq(£) satisfies V 2 Uq = where for brevity we have set a — 1, 
complemented by U$ = on the hole and periodic (anti-periodic) boundary condition 
on the edge of the cell. From the uniqueness properties of the Laplacian Uq = is 
the unique solution and hence usual homogenization cannot find the asymptotics at low 
frequency in this Dirichlet case. 



2.2. Repeated eigenvalues: linear asymptotics 

Repeated eigenvalues are commonplace in specific examples and, as we shall see, 
are particularly relevant to Dirac-like cones (Chan et al. 2012) that have practical 
significance. If we assume repeated eigenvalues of multiplicity p the general solution 
for the leading order problem becomes, 

u = f«\X)ulf\Z;n ) (19) 

where we sum over the repeated superscripts (I). Proceeding as before, we multiply 
equation (7) by Uq , subtract ui((aU^:)^ + VIqpUq ) then integrate over the cell to 
obtain, 

/o ( l / I ut\2aU& + a,M l) )dS + n?/^ / f g pU^U^dS = 0. (20) 

There is now a system of coupled partial differential equations for the /q and, provided 
Qi 7^ 0, the leading order behaviour of the dispersion curves near the Qq is now linear 
(these then form Dirac-like cones). These coupled partial differential equations on the 



long-scale now replace ( 15 ) near these frequencies 



For the perfect lattice and Bloch wave problem, we set / 



m 



obtain the following equations, 



Ki . „■»,-, \ i(l) 



1 ■"■jml ~r *L\"ml j Jo 



where, 



and 



i-jml 




U [m) (2aU { 



(0 



0, for m 



a,tM l) )dS 



1,2,..., 



P 



/o exp(iKjXj/e) and 



(21) 



B 



ml 




P U^ut ] dS 



The system of equations (21) is written simply as, 
CF = 0, 



(22) 



(23) 



(24) 
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Table 1. The coefficients an, necessary in equation (25) for the two Dirac-like cones 
at point r of the Brillouin zone for a doubly periodic array of constrained points. 

with Cu = £lfBu and C m ; = iKjAj m i/e for / ^ m. One must then solve for 
Qf = ±y/aijKiKj/e when the determinant of C vanishes and the asymptotic relation 
is, 

^ ~ ^0 ± yfry/OlijKiKj. (25) 

If Qi is zero, one must go to the next order and a slightly different analysis ensues. 

2.3. Repeated eigenvalues: Quadratic asymptotics 

Assuming that fii is zero, U\ = Jqx^xI ( we a S a i n sum over a U repeated (I) superscripts) 
and advance to second order using (|8J). Taking the difference between the product 
of equation (8 
elementary eel 



) with Uq and W2((aE/o,&),& + ^IpUo) and then integrating over the 
gives 



A% iXi J I aUtH^S + ti% kXj J f s ut\2aU^ + a^)dS 

+ nifS ] J J pUt>u!PdS = 0, for m=l,2,...,p (26) 

as a system of coupled PDEs. The above equation is presented more neatly as 

fo% iX Arni + fo% k xPkjrni + ^ l) B ml = 0, for m = 1,2, ...,p. (27) 

For the Bloch wave setting, using /q (X) = /q ' exp(iKjXj/e) we obtain the following 
system, 

^-^A ml -^-B kjml + n 2 2 B ml ^ff = 0, for m = l,2,...,p (28) 

and this determines the asymptotic dispersion curves. 

3. Effective dispersive media 

The major goal of HFH is to represent the periodic medium of finite extent using an 
effective homogeneous medium and the main result of this article is in achieving this, 
and for clarity this is summarized in this section back in the dimensional setting. 



Transforming equation (15) back to the original Xi = X{L coordinates and using 



the solution of f^ we obtain an effective medium equation for f that is, 

Tjo^x) + ^-^/o(x) = 0. (29) 
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Figure 4. The dispersion diagram for a doubly periodic array of square cells 
with circular inclusions, of radius 0.4, fixed at their inner boundaries shown for the 
irreducible Brillouin zone of Fig. [3] The dispersion curves are shown in solid lines and 
the asymptotic solutions from HFH are shown in dashed lines. 



The nature of equation (29 ) is not necessarily elliptic since Tn and T 2 2 can take different 
values and/or different signs. In the illustrations herein the Tij = for i ^ j. The 



hyperbolic behaviour of equation (29) yields asymptotic solutions that describe the 



endoscope effects where, if one of the Tu coefficients is zero as often happens near the 
point M of the Brillouin zone, waves propagate only in one direction. Note that when 
the cell has the adequate symmetries the dispersion relation near point G(0, tt/2) is 
identical to that near point X(n/2, 0) which explains the existence of two orthogonal 
directions of propagation instead of only one. 

At Dirac-like cones the linear behaviour of the effective medium yields an equation 



slightly different from (29). Regarding the zero radius holes, presented in the following 



section 4.3, it consists of a system of three coupled PDE's that uncouple to yield one 



?(0, 



identical PDE for all /q s that is of the form, 
/o^ < (x) + /3 ( " 2 "" g)2 /o(x) = 



P 



0. 



(30) 



where /3 is a coefficient equivalent of the 7L- but this time is the same for all combinations 
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0.698841854976085 
7.867675441589871 

0.3230 

4.9681 

4.775527931399718 

-4.279843054160115 



T22 
0.698841854976085 

7.867675441589871 

-8.998 

14.2899 

4.775527931399718 

-4.279843054160115 



Q 
1.700916617638699 
3.361627184739501 
3.632730109763024 
3.632730129014608 
4.148661549527329 
4.877003447953185 



Table 2. The first six standing wave frequencies for in-phase waves at T, cf. Fig. |1J 
together with associated values for T\\ and T22. Symmetry between T\\ and T22 is 
breaking when the multiplicity of the eigenvalue is greater than two. Negative group 



velocity is demonstrated by the negative sign of both X^ 's. 



of i and j. The quadratic mode that emerges in the middle of the linear ones follows 



equation (29) 



4. Dispersion curves 

We now generate dispersion curves for circular and square holes and verify the HFH 
theory versus these numerically generated curves; this is a good test of the approach 
as the dispersion curves contain changes in curvature, coalescing branches and modes 
that cross. These dispersion curves can then be used to interpret the physical optics 
phenomena seen later. Most vividly the dispersion curves also show the zero frequency, 
and other, stop bands. 

4-1. Circular inclusions in a square array 

The circular inclusions are arranged in a square array and each elementary cell is a square 
of side 2, we consider large holes initially and then how the curves change as the radius 
decreases. Ultimately we consider infinitesimal radii and perhaps remarkably obtain 
explicit solutions. In almost all other cases one has to use fully numerical techniques, 
or semi-analytical methods such as multipoles (Zolla et al. 2005, Movchan et al. 2002). 



4.1.1. Large circular holes The full dispersion curves are computed numerically using 
COMSOL (finite element code (COMSOL 2012)) and then asymptotically with the HFH 
using the equations developed in section[2J Fig. [4] illustrates the typical dispersion curves 
versus the asymptotics for wavenumbers, on the specified path of Fig. [3j for the example 
of a hole with radius 0.4. There is, as expected, a zero-frequency stop-band and the 
lowest branch is isolated with a full stop-band also lying above it. At some standing 
wave frequencies two modes share the same frequency, for instance the third and fourth 
modes at T; these modes are asymptotically quadratic in the local wavenumber. Table 
[2] shows the T U ,T 2 2 values for point T of the Brillouin zone. Note how T n = T 22 for all 
single eigenfrequencies, but that symmetry breaks for the double roots. Moreover the 
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c; 




Figure 5. The solid line dispersion curves are from FE computations (circle of radius 
0.4) and the dashed are from the HFH asymptotics. These curves become virtually 
indistinguishable when one takes larger and larger macro cells, as this amounts to 
perturbing away from an increasing number of standing waves on the diagrams, which 
is the essence of the folding technique. In the limit of an infinite macro cell one would 
perturb away from a dense set of points on the diagrams, and their reconstruction 
would become perfect (Cherednichenko 2012). 



signs of the Tn and T22 naturally inform one of the local curvature near V. Physically, 
this tells us the sign of group velocity of waves with small phase-shift across the unit 
cells, and thus whether or not they undergo backward propagation, which is one of the 
hallmarks of negative refraction. 



4-1.2. Folding technique for reconstructing the dispersion curves asymptotically An 
apparent deficiency of the approach we present is that it requires known standing wave 
frequencies, and associated eigenstates, at the band edges and that the asymptotics 
may be poor at frequencies far from these standing waves; until now HFH, applied 
to the dispersion curves, has been limited to obtaining asymptotic solutions near the 
standing wave frequencies (Craster et al. 2010). This uses an elementary cell containing 
a single hole and the irreducible Brillouin zone associated with it, however by using 
larger macro-cells containing four or more holes, and foldings of their resultant Brillouin 
zone, one can extract the asymptotics at other positions in wavenumber space thereby 
considerably enhancing the asymptotically coverage, see Fig. |5j As an example, we take 



an elementary square cell of length 2 with circular holes of radius 0.4, section 4.1 and 
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0.979719694474189 

1 
12.439613406598808 

-16.77 
18.856092144598160 



T22 
0.979724590049377 

-22.34 
12.439765169617527 

18.77 
18.855405067087230 



Q 
0.624563144241832 

7T 

3.388815608872719 
4.670779386112907 



Table 3. The first five frequencies for in-phase standing waves and their respective 
Tij coefficients for a doubly periodic array of square cells with zero radius holes (see 
Fig. [6|. Equal T\\ and T22 coefficients are for the single multiplicity eigenmodes. The 
second and fourth set of T^ are for the quadratic modes of the Dirac-like cones. 



Fig. § 

To obtain asymptotic solutions at the wavenumber positions, I(7r/4, 0), J(n/2, 7r/4), 
M'(7r/4, 7r/4) of Fig. gb), bisecting the three segments [TX], [MX] and [TM] of the 



Brillouin zone we can use the equations (15, 16, 17) again, section 2.2, but now with 
different elementary cells. 

Let the Brillouin zone for a square macro-cell composed of four holes be T'M'I. At 
r'(0, 0) periodic (in-phase) conditions on the macro cell boundaries yield all in-phase 
and out-of-phase modes of the elementary cell containing a single hole. This elementary 
cell's [TM] segment is folded on point M' . At M', where anti-periodic (out-of-phase) 
conditions on the macro cell are applied, the dispersion curves contain repeated roots, 
which are linear, and their slope is equal to the asymptote of the dispersion curves of 
the elementary cell setting at point M' . Proceeding in a similar manner one can fold the 
Brillouin zone at will by considering macro-cells composed of simple unfoldings of the 
elementary cell. Asymptotics can be obtained for wavenumber positions of the Brillouin 
zone equal to all ratios of the elementary cell's Brillouin zone wavenumber points. These 
ratios are equal to 1/n where n is related to the number of elementary cells of one hole 
contained in the macro cell in each dimension. 

We now illustrate this by obtaining asymptotics at points M', I and J which 
respectively correspond to the following macro cell settings, a square macro cell 
composed of four elementary cells with anti-periodic boundary conditions, a rectangular 
macro cell composed of two elementary cells in the £1 direction with anti-periodic 
boundary conditions on £1 but periodic on £ 2 and finally a rectangular macro cell 
composed of two elementary cells in the £2 direction with anti-periodic conditions in 
£1 and £ 2 . The asymptotics at points M', I, J as well as the usual points T, M and X 
are illustrated in Fig. [5j 



4-2. Circular inclusions of small radius 

As the hole radius decreases, the dispersion curves gradually transition toward those of 
the infinitesimal holes shown in Fig. [6] The zero frequency stop-band is generic, and a 
feature forced by the Dirichlet boundary condition, however, the lowest dispersion curve 
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Figure 6. The dispersion curves for a doubly periodic array of constrained points: 



Solid lines are from exact dispersion relation (32) and the dashed lines are HFH 



asymptotics. Additionally, we draw two sets of dotted lines emerging from T and 
X which are folded light lines. 



is no longer isolated and triple crossings created by repeated roots occur. 



For the larger hole radius, section |4.1[ all the modes are locally quadratic close 
to the standing wave frequencies. As the wavenumber moves away from those specific 
Brillouin zone points the behaviour of the dispersion curves becomes locally linear. The 
formation of triple, and more, crossings as the radius tends to zero, are created by 
the linear behaviour of the dispersion curves far away from points T, M and X moving 
toward those points to replace the quadratic behaviour of all the multiple modes, except 
one, as seen in Fig. [6j Triple crossings are illustrated at the second and fourth standing 
wave frequencies at point T and the first standing wave frequencies at point M, where 
the third frequency at point M yields a hep-tuple crossing. Reassuringly, the asymptotic 
theory captures these multiple modes whether the crossings are double, triple, or more 
as is illustrated in Figs. [4j [6] where the circular holes are respectively of radius and 
0.4. For the multiple crossing points with multiplicity higher than two, both linear and 



quadratic dispersion curves arise. One can proceed, as in section |2.2[ to obtain the 
different Qi coefficients and one of them is zero; proceeding to higher order and apply 



equations (26) to (28) gives the quadratic behaviour of the middle mode emerging from 
the triple crossings. Again the T ti contain the critical information about the local 
behaviour, some typical values along T are given in table |3j for the quadratic curves; 
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they are naturally identical for isolated modes but describe the inherent anisotropy for 
the other cases. One can then immediately see whether the material will behave with 
directional preferences at specific frequencies. 

Dispersion diagrams of zero radius and for 0.1 radius (not shown) geometries are 
almost indistinguishable by eye, but apparent triple crossings in the latter are actually 
a double crossing together with a single mode that has almost the same standing wave 
frequency. The separation of the two almost touching standing wave frequencies depends 
on the size of the inclusion and tends to zero as the inclusion shrinks down to a point: 
Fig. [7] illustrates a close-up of a triple point comparing inclusions of radius and 0.1. 
This is actually important for the Dirac-like cones and associated effects. 

4-3. Zero radius, Fourier series and Dirac-like cones 

A doubly periodic array of constrained points (circles of zero radius) is an attractive 
limit and has an immediate solution in terms of Fourier series as 

oo oo -i7rN-£ 

«&,&) = e te * £ £ 7 y^~ ( ^^ (31) 

ni=~oori2=—oo 

with k = (ki, k 2 ) an d N = (m, n 2 ) and Q and k are related via the dispersion relation 

oo oo 1 

D(k,Q)= V Y] -, - 9 — -. - 9 — 7^ = 0. (32) 

v ' ^ ^ (ki - 7mi) 2 + (k 2 - im 2 ) 2 - 9? v ' 

ni=—oo ri2=— oo v ' v ' 

This is valid provided the double sum does not have a singularity, which would occur 
whenever 

£l = a/(ki - nni) 2 + (k 2 - 7in 2 ) 2 . (33) 

These relations from the singularities also play a role in the dispersion picture; 
they correspond to the Bloch states of a perfect medium without any constraints 
(a homogeneous isotropic medium). However, in some circumstances they also, by 
serendipity, exactly satisfy the constraint u = at the centre of each cell and therefore, 
in those cases, are simultaneously dispersion curves; this is the origin of the degeneracy 
seen in the multiple crossings (occurring at Dirac-like points) of the dispersion curves 
of Fig. |6j It is also clear from the dispersion relation that Q = is not a solution at 
k, = and hence that there is a zero-frequency stopband and conventional long-wave 



homogenization is doomed to failure cf. section 2A_ 

A set of dispersion curves are shown in Fig. |6j The Bloch states for the perfect 
medium free of defects lead to light lines folded at the edges of the Brillouin zone and 
those emerging from T are pertinent to the current discussion. These light lines intersect 
at the multiple crossings, which are called generalized Dirac-like points: Generalized 
Dirac-like points occur for a given set of frequencies at all three edges of the Brillouin 
zone, namely for points T, M and X. It is useful to find criteria for their multiplicity: 
At T the generalized Dirac-like points occur for Q 2 = n 2 m such that m = n\ + n\ with 
m integer; the multiplicity depends on the ways in which n± and n 2 can form m. From 
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elementary number theory m can be equal to the sum of two squares in more than one 
way, i.e., m = 50 = 5 2 + 5 2 = 7 2 + l 2 or the sum of two squares can also be a perfect 
square, i.e., m = 25 = 5 2 = 4 2 + 3 2 . Equation u tXiXi + tt 2 {ti\ + n\)u = together with 
the appropriate boundary conditions describes the generalized Dirac-like points and an 
independent set of solutions is formed from different possibilities of sin(7r(njX + Tijy)), 
sm(n (niX — rijy)), cos(ir(niX+njy))-cos(ir(n p x — nkfj)) wherei,j,p,k G {1,2}, i ^ j and 
p y^ k. The multiplicity can be equal to 3, 7, 11 or more depending on m. For example 
Q 2 = 7r 2 with Hi = 1, rij = has multiplicity 3 where Q 2 = 57r 2 has multiplicity 7 and 
Q 2 = 50tt 2 has multiplicity 11. This is true for solutions at T and M but at X one 
can obtain a singularity that accepts only one eigensolution. Due to the non-symmetry 
of the general sum at X that renders the denominator of the Fourier series singular, 
m = (1/2 — rii) 2 + n 2 , it is possible to obtain solutions of multiplicity one as in the 
case of Qq = tt/2 where the only eigensolution respecting the boundary conditions is 
sin(7rx/2). 

Given the exact solution one can, in these special cases, generate the asymptotics 



by hand. For the asymptotics we construct the coefficients using U$ from (31) and 
deduce Ui = (Uu, C/i 2 ) as 

u» = 2*r* £ £ T i^ph ^^ (34 ) 

*— ' ^-^ («i — 7mir + («2 — 7TO2) — " 

ni=— 00 ri2=— 00 Lv y v ' J 

and thus one can extract the Tij in (15) by doing the integrals by hand, the off-diagonal 
terms are zero, and T n and T 22 are 

Eoo y^°° (ki— 7rni) 2 

ni=-oo 2^n 2 =-oo [( Kl _ 7rni )2 + ( K2 _ 7rn2 )2_Q2j3 

T n = 1 - 4 y^ ™ 1 (35) 

2—in\=— CO Z^rt2=— OO [( Kl _ 7rni )2_|_( K2 _ 7rn2 )2_Q2]2 

and an identical equation for T 22 but with 1 and 2 interchanged, where (k±, « 2 ) are (0, 0), 
(7r/2, 0) and (tt/2, tt/2) at the edges of the Brillouin zone at the respective points of T, 
X and M. These asymptotics only apply at the isolated, non-repeating roots, and are 
shown in Fig. [6j 



The linear cases at the triple crossings are easily extracted from section |2.2| For 
instance, for the repeated case at Q Q = tt there are three linearly independent solutions 
with 

u = f {1) (X) sin tt^ + f {2) (X) sin 7r& + / (3) (X) [cos ir& - cos tt^] (36) 

Following through the methodology one arrives at 2Vt\ = (2tt) 2 (k, 2 + k 2 ,) which gives 
the two linear asymptotics, and similarly at the other points. Table [l] summarizes the 
coefficients for the first two Dirac-like cones at point I\ The asymptotics are shown in 
Fig. [7] showing that the linear behaviour is perfectly captured. 

4-4- Square inclusions 

Our methodology is not limited to circular inclusions and is easily applied to any shape, 
we briefly consider square inclusions taking a square of side a/2 and see in Fig. |8| 
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Wavenumber 



Wavenumber 



Figure 7. Panel (a) shows a close up of the Dirac-like point and its asymptotics from 
figurepl The values for the asymptotic coefficients of equation (25 1 are an = c*22 = 27r 2 
and for equation (181 they are Tu — 1 and T22 = —22.34. Panel (b) shows the same 



region for a slightly larger hole radius of 0.1 where one can see the merging of the double 
mode with the mode near -n to form a triplet as the radius decreases to zero. HFH 
captures both cases as is seen from the asymptotic curves. As the radius decreases the 
top and bottom quadratics become steeper and ultimately reach a linear behaviour. 



G 




-1.5 

M 



0.5 1 1.5 

Wavenumber 



2 2.5 

x 



M 



Figure 8. The dispersion curves for square inclusions, of side \/2, from numerical 
simulation (solid) and from HFH theory (dashed). The HFH and numerics for the 
lowest two curves are virtually indistinguishable. The crosses on the frequency axis 
are frequencies predicted by solving Hclmholtz's equation in a waveguide consisting of 
a rectangle with Neumann data on one side and Dirichlet data on the other sides, see 
also Figure. 10. Their values in increasing order are 5.59, 6.22, 7.14 and 8.26. 
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Figure 9. The dispersion curves for a square of side v2 from numerical simulation 
for (a) a square rotated by 7r/8 and (b) by 7r/4. The crosses on the frequency axis in 
(a) are frequencies predicted by solving Hclmholtz's equation in a waveguide. Their 
values are 7.51, 8.24 and 9.33. 



that the dispersion curves are not unlike those of the large cylinder in Fig. [4[ an 
isolated lowest mode separating a zero frequency stop-band from a wide stop-band, yet 
another stop band arises near Q = 6. The asymptotics from HFH again reproduce the 
behaviour near the edges of the Brillouin zone and can be used to construct effective 
HFH equations. The HFH results in Fig. [4] are almost completely indistinguishable all 
along the lowest two branches and are highly accurate near the edges of the Brillouin 
zone for the other branches. Interestingly, the orientation of the square matters strongly 
and rotating it progressively flattens the lowest (lower frequency) dispersion curves 
until ultimately, for a rotation of 7r/4, they become completely flat leading to resonant 
states: The rotation gradually isolates each piece of material from its neighbours leaving 
resonant blocks that vibrate with the eigenfrequency of Dirichlet squares of side y/2 



n 2 



[7TITC 



2) 2 + (mr/y/2) 2 with m,n non-zero integers. 



4-4- 1- Standing wave modes It is noticable in Figs. [8] and |9^l that completely flat 
modes exist which represent standing waves for a whole range of Bloch wavenumbers. 
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1 





Figure 10. The eigenstates of the first four flat bands at the respective frequencies 
of 5.56, 6.11, 6.96 and 8.03. The left of every panel shows FEM computations of the 
eigenstate, as noted in the text a waveguide model can be developed to approximate 
the field and eigenvalue and these waveguide counterparts are shown alongside. 



The eigenfunctions, for the non-tilted case, for the first four standing wave frequencies, 



are shown in Fig. 10; it is clear that the field is concentrated along parallel and opposite 
sides and this suggests a simple equivalent waveguide model can be constructed. Taking 
a rectangle with a Neumann condition on one side and Dirichlet conditions on the 
other three sides where the Neumann condition on one side is necessary to cover the 
appropriate symmetry. The length and width of the rectangle depends on the tilt of the 
square hole. For the non-tilted case the length is taken to be the length of a cell and 
the width as the half width of the cell minus the inner square. The resulting frequency 
estimate Qe reads, 



n. 



IX 



In- 1 
/ — w 



+ 



for n,m GN* 



(37) 



where / and w/2 are respectively the effective length and width of the waveguide. For 
the non-tilted square hole lo = 2 and wo = 2 — a/2 where for a tilt of 7r/8 they reduce 
to the values of /„./« = = 1.6051 and w^/s = 0.4335. Equation (37) then leads to the 



<■ tt/8 



'x' crosses placed on Figs. [8] and [9^, which are good predictions of the standing wave 
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(b) 



(d) 







Figure 11. Dirac-like cone versus cloaking: (a) Perfect transmission through a column 
of seven constrained evenly spaced inclusions of radius 0.01 (array pitch 2) with periodic 
conditions on either sides and Perfectly Matched Layers (PML) on top and bottom, 
and a plane wave incident from above at frequency J7 = 2.215; (b) Nearly perfect 
transmission for a plane wave at frequency fio incident from above on a finite array 
of 140 constrained inclusions with PML on either sides of the computational domain; 
(c) Same as in panel (b) when a clamped rectangular obstacle is placed inside the 
array; (d) Plane wave incident from above on the clamped obstacle at frequency i7 
for comparison; The slightly reduced amplitude in forward scattering in (c), but the 
lack of backward scattering, is a hallmark of cloaking. 



frequencies, and the approximate eigensolutions are shown in Fig. 10 The precision of 



the estimates degenerates as the frequency increases as the field near the ends of the 
imaginary rectangle gradually leaks into the undisturbed ligaments. 

5. Lensing, guiding and cloaking effects 



We now connect the HFH theory with the exciting phenomena that are topical in 
photonics. 
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Figure 12. Panel (a) shows nearly perfect transmission for a plane wave at frequency 
S7 incident from above on an effective medium with properties computed by HFH 
for a frequency near fl = 2.215 with PML on either sides of the computational 
domain; A Plane wave incident from above on the clamped obstacle at frequency 
S7 f° r comparison. 
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Figure 13. Lensing effects in a PC (tilted array) of small Dirichlet holes: (a) Lensing 



for a line source inside the PC at frequency ft — 2.7; (b) Using equation (29) the PC 
is replaced by an effective medium in order to yield the same effect as in (a). 



5.1. Cloaking effects near Dirac-like cones 

As noted in recent articles (Chan et al. 2012) a curious phenomenon occurs in photonics 
near Dirac-like points which are the triple crossings shown in, for instance, Fig. [7| Let 
us consider the lowest frequency triple point, near wavenumber M, shown in Fig. [6j for 
the full finite element numerical simulations we actually consider very small cylinders of 
radius 0.01 in a square array. The triple crossings occur because the dispersion relation 
of a perfect medium (no cylinders) consists of light lines folded in space (one folded light 
line is shown in Fig. [6]). The folded light lines happen to satisfy the boundary condition 
for the infinitesimal point holes, and are perturbed slightly for finite cylinders (as in Fig. 
[7| ; the middle curve passing between the Dirac-like cone is created by the inclusion and 
so these effects co-exist. Importantly, this means that there is near perfect transmission 
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Figure 14. Guiding and lensing effects in a PC (tilted array) of small Dirichlet holes: 
(a) Omni-directive antenna for a line source inside the PC at frequency fl = 7r/2; (b) 



The PC of (a) is replaced by an effective medium obtained by HFH with equation (29) 
in order to obtain the same physical effect; (c) Endoscope for a line source at £1 = 2.7; 



(d) The PC of (c) is replaced by an effective medium using equation ( 30 1 to obtain the 
same physical effect. 



through an array of cylinders close to these triple crossing frequencies as shown in Fig. 
IIVa,b) as the incoming plane wave, at f2 = 2.215, does not see the inclusions. Inserting 
a large Dirichlet hole within the array leads to virtually no reflection, the plane wave 
is hardly reflected, and its transmission is nearly perfect Fig. [lite) except for a slight 
shadow zone in clear contrast to the uncloaked Dirichlet hole Fig. [TlTd); If we remove 
the array of small holes, the defect is strongly scattered by the wave. 



Replacing the array of holes with HFH is shown in Fig. 12 the quasiperiodic 
medium is replaced by a homogeneous medium with the effective properties given by 
HFH. The excitation is at a frequency Q = 2.12, close to the standing wave frequency 
Qo. Inserting the values for I = 1, Q , Q and j3 = 2/7T 2 into equation (30) we obtain the 
effective continuum equation, 

/o,^(x) + 0.0393/ (x) = 0. (38) 

Simulations show qualitatively the same effects, one has transmission through the slab 
with plane waves incoming. 



26 




I 




Figure 15. Lensing effects in a PC (tilted array) of small Dirichlet holes acting as a 
flat lens near a frequency of tt/2: (a) A line source near the top of the PC yields an 
image on the other side. The PC is composed of 96 very small holes arranged in a 
tilted square array with an angle of 7r/4. (b) shows clearly the same effect happening 
when the PC is replaced by a anisotropic effective medium obtained through HFH. 



5.2. Lensing and wave guiding effects 

Returning to the introduction we show in Fig. J2J(a,c) a PC composed of an array of 196 
holes with radius 0.4; the corresponding dispersion curves are shown in Fig. |4| 

We now interpret these strongly anisotropic beams and how the HFH represents 
this. In Fig[2](a,b) we choose frequency Vt = 1.98, which from Fig. Elis near the standing 
wave frequency Qq = 1.966 at point X(tt/2, 0). The reciprocal space is only shown for a 
portion of the irreducible Brillouin zone and one can use reflections of the triangle chosen 
to fill the entire square in the Brillouin zone; due to these internal symmetries the effect 
is not only due to the first mode of Fig. El at point X(0, 7r/2) of the reciprocal space, 
but also to the first mode, at point G(tt/2,0) (not shown). The effect is reproduced in 
Fig. [2] (b) by combining two effective medium equations, one for point X and one for 
point G, each one is responsible for a single diagonal, and the HFH equations are 



-1.4778/U^ (x) + 0.8837/0,^ (x) + (fi 2 - fig)/ (x) = 0, 
0.8837/ 0i:Ci:l , (x) - 1.4778/0,^ (x) + (fi 2 - fig)/ (x) = 0, 



(39) 
(40) 



with (Q 2 — Qq) = 0.0552. Note that these effective equations have opposite signs in the 
coefficient, and so the governing equations are hyperbolic and not elliptic, the result 
is that these direct energy along rays or characteristics and the angle is given by the 
relative magnitude of the T\\ and T22; in this case roughly equal with the energy directly 
along the diagonals. 

To illustrate this further in panel (c) of Fig. [2] we show a St George's cross effect 
obtained near the standing wave frequency Qq = 2.744. By homogenizing the PC and 
exciting it at a frequency of Q = 2.75 we obtain two effective medium equations, one 
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for X and one for G that are, 

3. KM/o,^ (x) + 0.085/0^ (x) + (tt 2 - flg)/o(x) = 0, (41) 

0.085/0,^^ (x) + 3.1094/0,^^ (x) + (fi 2 - ng)/ (x) = 0. (42) 

Each equation is responsible for a single line of the cross, the orientation of this cross is 
explained by the almost zero coefficient, in each effective equation, that does not permit 



propagation in the vertical and horizontal direction for the respective equations (41) 



and (42). The coefficient in front of /o( x ) is equal to 0.033. It is clear therefore that 
the strong anisotropy that is witnessed in full FE numerical simulations has its origins 
in the size and sign of the Ty coefficients from HFH and the effective equations can be 
used to gain quantitative understanding and predictions. 

Another striking application involving wave propagation and an array of small holes 



of radius 0.01 is shown in Figs. [13j [14] , [15] wherein we have tilted the array through 
an angle 7r/4. The dispersion curves in Fig. M are for holes of radius exactly zero, but 
provide very good comparison for this small holes case. 

We obtain a highly-directed emission at frequency Q = 1.6 near the standing wave 



frequency f2o = 7T /2, shown in Fig. |I4[a) with its effective medium counterpart in 14 b); 
in the dispersion curves this is the lowest standing wave frequency at X. There are two 
effective medium equations at this frequency which yield, 



11.42/ ,* 1 * 1 (x) + /o,Wx) + (ft 2 - 


-fi 2 )/o(x) = 0, 


(43) 


fo,xixi(*) ~ H.42/ 0)XaXa (x) + (tt 2 - 


-fi 2 )/o(x) = 0, 


(44) 



where the difference of the frequency squares is equal to 0.0926. For forcing within 
the slab one again has a strongly anisotropic response, the waves within the medium 
forming a cross-shape which, when it strikes the edge of slab all gets radiated out into 
the surrounding medium. Placing a source outside the slab, again at a frequency near 
tto = 7r/2, Fig. [TB^a) shows a PC that behaves like a flat lens. In panel (b) of that Fig. 
the effective medium yields the same effect and is also governed by two equations of the 



form of (29) with ly coefficients equal to the ones of equation (44). The frequency at 
which the lensing effect is obtained is Q = 1.56 and so the coefficients in front of /o( x ) 
are equal to —0.0338 in this case. 

Figs. [j~3~|(a) and U4Tc) show lensing effects where the PC is excited at fi = 2.7, 
this frequency is chosen as follows: When the array is rotated the relevant light line in 
the surrounding material, and its folding, emerge from the point M and are shown in 
Fig. M we see that the energy from the source couples (as the light line crosses the full 
dispersion curves at Q = 2.7) into one of the linear Dirac-like cone modes that passes 
through Q = it; this is the standing wave frequency of the Dirac-like cone at I\ We 
therefore model the effective medium with HFH using the analysis near Q = tt, and it 



yields very similar results in Figs. 13 b) and (141(d), by homogenizing the PC with an 
effective material governed by equation, 

/o,^(x) + 0.337l/ (x) = 0, (45) 
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Figure 16. Lensing and guiding effect through a PC composed of square Dirichlct 
holes which we vary their tilt angle: (a) Lensing effects in a PC composed of Dirichlct 
square holes tilted by 7r/8 near the respective frequencies of 7.7 and 8.25 for the top 
and bottom figures (Fig[9^). (b) The dispersion curves in Fig. [9b show the band 
structure is composed of many band gaps separated by standing wave modes, so that 
the structure acts as a perfect reflector for any frequency in order to guide light, (c) 
shows a perfect reflector for a source exciting the medium at frequency 4.95, located 
in the second band gap of Fig. [8J (d) shows the same effect as in (c) for a frequency 
of 5.0 although in the present the filled parts of material are excited but no energy 
escapes the PC and acts as an energy trap. 



obtained by simply replacing the source and standing wave frequencies together with 
f3 = 1/(2tt 2 ) in equation (g. 



Finally, Fig. |16[a) shows a unidirective PC for the first two flat modes of the band 
diagram in Fig. [9Va). The standing wave frequencies of the two flat modes in question 
are Qq = 7.6 and Qo — 8.28. The sources are excited near these frequencies at Q — 7.7 
and Q = 8.25. HFH yields effective medium equations for each of the modes that are 
respectively, 



9.3649/c,^ (x) + (Q 2 - fig)/ (x) = 0, 



(46) 
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and 

-45.8434/0,^ (x) + (Q 2 - ^)/ (x) = 0, (47) 



Equations (46) and (47) clearly show the unidirectivity of the two effective media. 
Similar equations with interchanged coeffcients stand, for the symmetric location of 
the Brillouin zone G(0,7r/2). Panels (b), (c) and (d) contain PCs with rotated square 
holes by an angle of 7r/4 with an unsual band diagram seen in Fig. [9tb). Panel (b) 
shows the guiding of a wave by means of perfect reflection for a line source excitated 
just outside the standing wave frequency of Q = 4.95 while in panel (c) the same PC but 
this time rotated by an angle of n/2 acts as a perfect reflector. Panel (d) finally excites 
the medium right on the standing wave frequency of Q = 5 and standing waves appear 
between the holes. Small gaps between the hole's edges make it possible for energy to 
pass from one cell to the other. 

6. Concluding remarks 

In this paper, we have presented a range of applications of the high frequency 
homogenization theory, in the context of zero-frequency stop band structures, which 
were previously thought to be non-homogenizable: HFH has succeeded in capturing 
the finer details of both dilute and densely packed photonic and phononic crystals. 
Striking physical effects such as cloaking via Dirac-like cones and directive antennas and 
endoscope effects have been provided, and fully explained, via HFH; importantly this 
is all shown to be related to the potentially anisotropic T^ coefficients that encode the 
local behaviour or their equivalent terms for multiple crossings. The range of unsolved 
homogenization problems in the wave community is vast, and this new method now 
opens the door to efficient simulation of multiscale structures. 

The current theory is limited to periodic, or nearly periodic media, however one 
can envisage extensions in a variety of different directions: In this paper we deal 
with photonic and phononic periodic structures, but it should be possible to adapt 
our results to quasi-crystals using a generalized form of the Floquet-Bloch theorem 
in upper dimensional spaces. Stochastic cases where the hole positions involve some 
random perturbation, say, might prove more arduous from a mathematical standpoint, 
although a physicist might argue the period would be replaced by the mean distance 
between the inclusions. The key point in order to relax the constraint of classical 
homogenization is the assumption that the Floquet-Bloch theorem is applicable at least 
to leading order. Whenever this can be done, or an extension of the Floquet-Bloch 
theorem can be used, HFH is applicable. Also of interest are diffusion problems in 
composites, classic homogenization is here again constrained by low frequencies (e.g. of 
a periodic heat source), whereas much of the exciting physics lies beyond the quasi-static 
regime again we hope that the insight provided by HFH will lead to progress in these 
directions. 
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